Glassy dynamics and domains: exact results for the East model 
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A general matrix-based scheme for analyzing the long-time dynamics in kinetically constrained 
models such as the East model is presented. The treatment developed here is motivated by the 
expectation that slowly-relaxing spin domains of arbitrary size govern the highly cooperative events 
that lead to spin relaxation at long times. To account for the role of large spin domains in the 
dynamics, a complete basis expressed in terms of domains of all sizes is introduced. It is first 
demonstrated that accounting for single domains of all possible sizes leads to a simple analytical 
result for the two-time single-spin correlation function in the East model that is in excellent quan- 
titative agreement with simulation data for equilibrium spin up density values c > 0.6. It is then 
shown that including also two neighboring domains leads to a closed expression that describes the 
slow relaxation of the system down to c ~ 0.3. Ingredients of generalizing the method to lower 
values of c are also provided, as well as to other models. The main advantage of this approach 
is that it gives explicit analytical results and that it requires neither an arbitrary closure for the 
memory kernel nor the construction of an irreducible memory kernel. It also allows one to calculate 
quantities that measure heterogeneity in the same framework, as is illustrated on the neighbor-pair 
correlation function and the distribution of relaxation times. 

PACS numbers: 64.70.Pf, 61.20.Lc, 52.35.Mw, 02.50.Ey 



I. INTRODUCTION 

Despite much progress in recent years, many aspects 
of structural glasses and undercooled liquids still escape 
a complete understanding^ 2 *^. Rather than studying 
the behavior of molecular glasses, one often investigates 
the behavior of simple models in the hope to capture 
the basic physics of such systems. The so-called East 
model is one of these simple models 5 . It is a classic Ising 
model with a trivial Hamiltonian in which the stochastic 
dynamics governing the change of spin leads to a com- 
plicated and highly-cooperative evolution of the system. 
In the East model, any spin has finite probability to flip 
from up to down or vice- versa only if the spin to the east 
of it (i.e., of higher lattice index) is up. Models of this 
kind are generally called kinetically constrained models 
or facilitating spin modelsSiSsSiiS. 

Such models are designed to mimioiMLiSii 3 . the kind 
of dynamics that take place in glasses^. Although 
the East model itself does not have a glass transi- 
tion at any finite spin density 5 , the decay of the sin- 
gle spin time-correlation function at low densities is a 
stretched exponentia&i^, following a functional form 
similar to that of the dynamic structure factor in glassy 
systemsi 5 ^*^. In fact, the typical spin relaxation time 
has been shown to behave as log(r) ~ log 2 (l/c), where 
c is the equilibrium density of up-spins, heuristically 
by Sollich and Evans^ and rigorously by Aldous and 
Diaconisi^, indicating an extreme slowing down for small 
c, suggestive of a transition at c = 0. More recently, the 



East model has been analyzed to examine the nature of 
dynamic heterogeneitiesiiA^ in frustrated systems. 

The typical relaxation times present in systems 
exhibiting frustration can be retrieved from time 
correlation functions. In glasses, mode coupling 
theory is one of the predominant descriptions for 
these correlation functions. Several approaches to 
mode coupling theories exist. In the context of 
glasses, that of Gotze and co-workers has been widely 

use( jl6| 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30 

The approach of Gotze and co-workers expresses the 
time correlation functions in terms of a memory ker- 
nel and then uses a certain ansatz in which the memory 
function is written in terms of the correlation functions 
themselves, yielding self-consistent equations. Oppen- 
heim and co-workers have addressed the formal points 
and justification of mode coupling theories in Fourier 
space^^ 2 ^^. Along similar lines, AndersenS has for- 
mulated a phase-space mode coupling theory for general 
fluids that leads to self-consistent equations for time de- 
pendent correlation functions. 

Mode coupling theory can be applied for both deter- 
ministic and stochastic systems. For some systems, such 
as the East model, the application of the most commonly 
used mode-coupling ansatz necessary to close the result- 
ing equation of motion for the spin autocorrelation func- 
tion leads to a spurious transition from an ergodic to non- 
ergodic phase at finite values of the spin concentration c, 
a result that is clearly at odds with simulation results. 
To analyze the failure of the closure approximation in 
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mode coupling theory, Pitts et al. have presented a di- 
agrammatic treatment that yields similar equations to 
those of mode coupling theories of the glass transition^. 
This treatment lends itself to a closure assumption that 
is very similar kind to the mode coupling theory of Gotze 
et al. Pitts et al. propose alternative closure assumptions 
to this mode coupling theory by summing subsets of dia- 
grams. The resulting predictions are in good agreement 
with simulations for high concentration of up-spins, but 
still decay too rapidly at lower spin density. Recent im- 
provements on this scheme have been carried out by Wu 
and Caoi^ based on a combination of matrix methods 
and mode coupling closure assumptions. 

Regardless of the precise formalism, mode coupling 
theories aim to describe slow long time behavior, and 
so should include all the "slow modes" of the system. 
For hydrodynamics of a simple fluid at moderate densi- 
ties, these slow modes are well-known, namely the den- 
sity, momentum and energy modes at large wave lengths. 
Correspondingly, when applied to Fourier modes of these 
hydrodynamical fields, mode coupling theory (e.g. in the 
formulation of Oppenheim and co-workers) yields well 
defined perturbative results when the correlation length 
is finite and the thermodynamic limit is taken. However, 
in the case of the East model, the relevant slow modes 
are less obvious, and previously developed mode coupling 
theories for this model^2iiii2& may have missed some of 
these slow modes. In fact the absence of some slow modes 
provides a possible explanation for the problem these the- 
ories have with describing the long time behavior at small 
c. 

The main purpose of this article is to identify the rel- 
evant slow modes in frustrated spin systems and to de- 
scribe the impact of the coupling of these modes to a 
specific spin variable. It is demonstrated that the ex- 
istence of slowly-relaxing spin domains of arbitrary size 
suggests a natural basis of slow modes in which quanti- 
tatively accurate but simple approximation schemes are 
easily formulated for many quantities of interest in the 
study of slow heterogeneous relaxation. 



II. THE EAST MODEL 

The East model^ is a linear chain of N spins, which 
are numbered from to N — 1, with each spin allowed to 
assume one of two values at any given time, here taken to 
be up or down. Occupation numbers rii are defined such 
that rii = 1 when spin i is up and if it is down. The 
static properties of this model follow from the Hamilto- 
nian 



H 



N-l 

H ^2 rii. 

i=0 



(1) 



Using the canonical distribution p ~ exp[— /3ff], the av- 
erage occupation per site if the system is at equilibrium 



at an inverse temperature /3, is found to be 

c=l/(l + e^). (2) 

As ftfi has little physical significance in the current con- 
text, the density c will be used as a parameter. If n 
denotes a spin state (i.e., a configuration of the N spins), 
the canonical equilibrium distribution can be written as 



P {n) = Y[ [cm + (i - c) (l - m)] 



(3) 



i=0 



where Eqs. Q and (J2J) were used, as well as the fact 
that rii is either zero or one. Note that each spin i has 
a probability c to be up (rij = 1) and 1 — c to be down 
{rii = 0). 

For the dynamics, consider the conditional probability 
density C/t(n, n') to be in state n at time t given that the 
system was in state n' at time 0, which satisfies^! 



£7 t (n,n') = 5^[w(n,n)l7 t (n,n / )-W(n,n)l7 t (n > n / 

n 

= ^£(n,n)[/ t (n,n'), 



(4) 



with C/o(n, n') = <5 nn '. Here, W(n,n')dt is the probabil- 
ity to make a transition from n' to n in a time dt. By 
definition, W(n, n) = 0. Defining 



A(n,t) = Y,A(n')U t (n',n), 

11 ' 

from Eq. it follows that 

i(n,i)=^£(n,n')A(n',t). 

or 

A(t) = CA(t) 



(•5) 



(0) 



where the Liouville operator £ is a linear operator on 
the 2 N dimensional Hilbert space of functions of n. The 
formal solution of Eq. © is A(t) = e ct A(0). 

For the East model, W(n, n') can be written as a sum 
over possible moves, i.e., spin flips of individual sites i: 

N-l 



W(n,n') = Wi(n,n'). 



i=0 

A flip of spin rii is possible only if rii+i = 1 (using the 
boundary condition that njy = 1), as the expression for 
Wi bears ou^: 

Wi(n, n') = [cS ni iS n '. + (1 - c)5 ni0 8 n <.i\ n i+ i JJ^nJ • 

Correspondingly, the Liouville operator can be written 

as 

N-l 

C(n, n') = [(! - c )^i (fyo " ( 7 ) 

i=0 

rij n . • 

3& 
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The equilibrium distribution in Eq. JSJ will serve as 
a weight for the inner product on the Hilbert space, 
i.e., the inner product of A(n) and B(n) is (A\B) = 
^] n p(n)A(n)B(n) = (AB). Only real quantities will be 
used in this paper, so there is no need to define a complex 
inner product (although this is straightforward). The 
time correlation function of A and B can now be written 
as (A(t)\B). When W(n, n')p(n') = W(ri ,n)p(n), i.e., 
when detailed balance holds, as it does in this model, C is 
Hermitian with respect to the inner product. In contrast 
to stochastic systems, in deterministic systems the Liou- 
ville operator is anti-Hcrmitian. It should be noted that 
the condition of detailed balance also guarantees that 
the limiting stationary distribution of the Markovian dy- 
namics is the equilibrium distribution J3"|. provided the 
underlying Markov process is ergodio^l. 

As remarked by Pitts et al£, different sites are not only 
statically independent [cf. Eq. Q], but also dynamically. 
To see this, consider the normalized single site fluctuation 



hi(t) 



i(t) 



^/cjl^c) 



(8) 



which satisfies (hf(t)) = 1 and (fii(t)) — 0. Note that 



Cfii = -m+ifti 

= —chi + \Jc(l 



c) rnn i+ i, 



(9) 
(10) 



using Eqs. J7J and ©. Thus, the time derivative of hi 
depends on the product of hi and hi + \. The deriva- 
tive of that product will in turn depend on hi+2, and 
so on. Thus, ?v(i) = e £t ?v involves only hi with 
i > i'. From the static independence of any site i 
and any other site i" with i" < i' < i, it follows that 
(hi»hi'(t)) — (hin)(hi>(t)) — 0. Because £ is Hermitian, 
also {hiu(t)hi') — 0. So all time correlation functions 
between different sites i' ^ i" are zero. 

Given the dynamical independence of different sites, we 
are interested in the nontrivial time correlation function 



C(t) = (h^h^O)). 



(11) 



In the limit N — > oo with i fixed, this is independent 
of i due to translation invariance. This single spin time 
correlation function in the thermodynamic limit is the 
main quantity of interest. 



III. PROJECTION OPERATOR TECHNIQUES 

In a mode coupling framework, dynamical equations 
are derived for the time correlation functions of slow 
modes in the system, which involves a memory ker- 
nel that is again expressed in terms of the correlation 
functions^S^S. The starting point is often a projection 
operator formalism. Here, a general setup will be pre- 
sented, which uses the Mori-Zwanzig projection operator 
formalism224£, to be specialized later. 



Let Ak be the slow modes of the system determined 
from physical arguments, with fc an index running over 
the slow modes. It will be assumed that (Ak) = (as 
could be achieved by subtracting the average), and that 
Ak are orthonormal, i.e., (Ak\A q ) = 5k q (as could be 
achieved by a Gramm-Schmidt procedure). For brevity, 
the Ak's are taken together in a vector A. In the projec- 
tion operator formalism, the component along A of any 
other physical quantity B is found using the projection 
operator 



VB = (B\A)-A = J2{B\Ak)A k , 



(12) 



where • denotes a vector product, i.e., a sum over fc, as 
indicated. Using the operator identity 



e ct _ e (l-V)Ct 



• f e CT VLe^ c ^dT, (13) 
Jo 



and Eqs. JHJ) and (|12fl . one can derive that 



A(t) = • A(t) + / U u (t - t) ■ A(r)dT + tp(t), (14) 



where ip(t) = e {1 ~ v)ct {l - P)CA and 



M E = {A\C\A} 
\ D (t) = (<p(t)\<p). 



(15) 
(16) 



Note that and M D (t) are matrices whose dimensions 
are equal to the number of slow modes and that M E 
contains only static information, while the memory ker- 
nel M D (t) involves the time correlation of the fluctuating 
force (f(t). 

Taking the inner product with A, Eq. 114fl yields an 
equation for the correlation function G(t) = (A(t)\A): 



G(t) = M E ■ G(t) + / M D (t - r) • G(r)dr. 
Jo 



(17) 



To solve this equation, a Laplace transform will be used, 
defined as 



G(z) = / e~ zt G{t)dt. 
Jo 



(18) 



Here and the following, we adopt the convention that 
quantities with a tilde (~) are z-dependent. The solution 
in Laplace space of Eq. (fTTfl is 



G = (zt — M E — M D ) 



D\-l 



where 



-ztf,i\D 



(t) dt. 



(19) 



(20) 
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IV. PHYSICS OF THE SLOW DYNAMICS 

If A in the previous section contained all the slow be- 
havior, then M D (t) would be a quickly-decaying function 
that could be replaced by a delta function in time and 
integrated over in Eq. 11171^ -^. Unfortunately, this is 
typically not the case because the projection operator 
(1 — V) only removes part of the dependence on A. For 
instance in a fluid, the long wave length modes of density, 
momentum and energy are slow because they correspond 
to densities of conserved quantities, but only at low den- 
sities is it enough to consider only these modes as slow. 
Extending the set A by multi-linear modes^^i, can help, 
and can be used to setup self-consistent equations which 
are exact in the thermodynamic limit provided there is a 
finite dynamical correlation length2L22iSiM. 

As is known from the extensive work of Gotze and co- 
workera 1 fi i 2n i 21 ,22,23,24,25,26,27,28,29,30 in deterministic sys- 
terns, such self-consistent equations can give rise to a 
glass transition. 

If, on the other hand, A is a complete set, then 
1 - V = and consequently ip{t) = and M D (t) = 0. 
In this case, the above formalism corresponds to writing 
Eq. 10 in a particular basis. This formulation is often ap- 
plied to the East model^iia. When working with a com- 
plete basis set, the set still has to be truncated at some 
level in practice. This introduces truncations errors, or, 
viewed alternatively, a nonzero M D (t). To get beyond the 
truncation problem, one makes an ansatz for the memory 
kernel in terms of the time correlation function of interest 
[here C(t)}, yielding a self-consistent equation. However, 
in stochastic systems, a glass transition will not be found 
if such an ansatz is used for the memory kernel, due to the 
Hermitian nature of £2*^. Rather, an ansatz needs to be 
used for the so-called "irreducible" memory kernel 36,42 . 
Then, a glass transition can be found for finite c in the 
East mode&S. However, simulations make it clear that 
there is no transition to a non-ergodic phase at a non- 
zero value of c. Somewhat better schemes to improve 
the ansatz have been developed sincei 4 *, but generally, 
they lead either to a transition or to a time correlation 
function that decays too quickly. 

Essential for the success of mode coupling theories for 
fluids at lower densities is the finiteness of the dynamical 
correlation length, which gives a cutoff length and cer- 
tain exact factorization properties in the thermodynamic 
Iimi1j2ii2£. However, for the East model, no such length 
scale exists, as is seen when one contrasts the result of 
the scaling of the decay tim o 18 i 19 vs. the diagrammatic 
approach of Pitts and Anderson^. When diagrams are 
truncated at a certain level, corresponding to taking into 
account only spins within a certain distance I, the mem- 
ory kernel becomes a polynomial in c of which highest 
power is c l , which means the typical time can scale at 
the slowest as c~ l . However, for c — ► 0, the timescale di- 
verges faster than any inverse power of aiS^. So clearly 
the spins at all positions, arbitrarily far away, need to be 
taken into account. 



A second way to see that the dynamic correlation 
length is unbounded is to note that while the static corre- 
lation length is zero, dynamically, the decay of the time 
correlation function C(t) is influenced by other spins. 
For example, if there is a large domain of down spins to 
the east of a given spin, that particular spin requires a 
long time to flip since all down spins in the domain must 
flip at least once. Hence, the decay is correlated with the 
existence of this domain, and the dynamic correlation 
length is therefore at least of the order of the size of this 
domain. But domains of all sizes exist and larger sized 
domains will contribute to the behavior of the time cor- 
relation function C(t) at longer times. Even if one is only 
interested in the bulk of the behavior of the time correla- 
tion function C(t), for which the relevant domains are of 
typical size 1/c, this size diverges as c — > 0. Therefore, it 
is no surprise that fixed spatial truncations do not work 
below a certain value of c and that mode coupling theo- 
ries using such truncations have problems to describe the 
long time behavior of the time correlation function C(t). 
For a different formulation of the importance of domains, 
see Garrahan and ChandleniiiiS and Wu and C&o—. 

Thus, physically, the origin of the slowness of the dy- 
namics seems to be related to the absence of a finite dy- 
namical correlation length and the existence of arbitrarily 
large domains of down spins. 

V. THE DOMAIN BASIS 
A. Single domains 

Consider the leftmost spin no in a semi-infinite chain 
of spins. East of this leftmost spin (i.e., at sites i > 
0), a domain of typical size 1/c filled with down spins 
exists. In the previous section it was argued that the 
presence of these domains is essential to the dynamics, 
so they should somehow be included. This is achieved by 
defining the (single) domain basis, which is composed of 
the orthonormal basis vectors 

Qo = ^(no-c) (21) 

and 

Qi(0) = ^^y(^o - c)( ni - c) (22a) 

= ■F7TT( n °-c)(l-n 1 )(n 2 -c) ( 22b ) 
Al(1J 

and in general 

k 

Qi(k) = -^Tjr (no ~ c) JJ(1 - nj )(n k+1 - c). (22c) 

where the normalization constants are to be chosen such 
that {Ql) = ([Qi(fc)] 2 ) = 1. Note that in Eqs. each 
factor (1 — rij) only yields a contribution when rij = 0, 



5 



k 



Qo Qi(0) 



TABLE I: Diagrammatic representation of the domain basis. 



i.e., when spin j is down. Thus, a consecutive sequence 
of such factors represents a down-spin domain. 

Comparing the basis vector in Eq. 1|21|) to ho in Eq. JHJ , 
which is normalized, it is seen that Qo should be equal 
to no, i.e., 



Zq = y/c(l - c). 



(23) 



Using a convenient diagrammatic representation we will 
later derive that 



Zi(fc) = c(l-c) 1+fe / 2 . 



(24) 



Other useful quantities will be the unnormalizcd versions 
of these basis vectors 



Qo 



(n - c) 



(25) 



h(k) = (no-c)JJ(l-nj)(nfc+i-c). (26) 

3=1 



Note that 



Qo 
Qi(fc) 



Qo/Zo 

Q 1 (k)/Z 1 (k). 



(27a) 
(27b) 



It is convenient at this point to introduce a diagram- 
matic notation for the unnormalized quantities Qo and 
Qi, depicted in Tabled In such diagrams, the horizontal 
direction represents the lattice. An open circle (o) at a 
certain position i denotes (n,i — c), a horizontal line (- ) 
denotes (1 — rii) and a product of consecutive (1 — nj)'s is 
represented as a longer horizontal line with the number 
of factors written on top of the line. When these dia- 
grammatic elements are touching, this denotes that they 
are on neighboring lattice sites. The first (i.e., leftmost) 
symbol in a diagram will always refer to site 0. Note that 
after a diagram for Qo and Qi has been evaluated, one 
still has to divide by the appropriate normalization to get 
the quantities corresponding to Q and Qi, cf. Eqs. (|2TJl . 

As an example, this diagrammatic notation will be 
used to show that Eq. (|24|) is the correct choice for Z\ (k) 
such that (Qi{k')\Qi(k)) is equal to Skk'- The product 
of Qi(fc) and Qi(fc') will be represented as the two dia- 
grams of Qi(fc) and Qi(fc') from Table [I] on top of each 
other: 



(Qi(fc')IQi(fc)) = 



(o— -o 



(28) 



The vertical displacement of the two parts in the diagram 
on the right hand side serves to distinguish one part from 
the other and to indicate that the diagram needs to be 
averaged. Next, note that since (o) = (n^ — c) = 0, any 



open circle at site i in the diagram has to overlap with an 
open circle or a horizontal line at the same site i from the 
other part in order to give a non-vanishing contribution 
to the average. In the diagram in Eq. (|28|l . this can only 
happen if k = k'. Thus, the domain basis is orthogonal. 
The case k = k' is represented by 



<Qi(fc)IQi(*)> = 8=^ =c 2 (l-c) 



k+2 



(29) 



In the last equality, we used that different parts of a 
diagram have different contributions depending on which 
diagrammatic elements on top align with which elements 
on the bottom, as shown in Table [H] (the diagrammatic 
element "•" will be introduced later). The contributions 
of the various parts are simply multiplied because each 
site is (statically) independent. From Eqs. (|27b() and (|29|) 
it follows that 



(Qi(fc)IQi(fc)) 



(Qi(fc)IQi(fc)) c 2 (i-c) 



k+2 



[Zi{k)f 



(30) 



With the definition of Zi(k), Eq. I|24|l . we immediately 
see that the Qi(fc) are properly normalized. In a similar 
way, one can deduce that (QolQi(fc)) = 0, thus establish- 
ing that the basis is really orthonormal. 

Taking the collection {Qo, Qi} for A in the projection 
formalism of Sec. IIIII the matrix M E in Eq. I|15l) which 
determines the dynamics, becomes 



(Qo|£|Qo> 
(Qi|£|Qo> 



(Qi|£|Qi> 



(31) 



where Qi without a value of k denotes the column vector 
(in the ket) or row vector (in the bra) composed of all 
Qi(fc). From Eqs. Q and @ it follows that 



CQ 
£Qi(0) 

CQ x {k> I) 



— (no — c)ni = — o» (32a) 

— (n — c)[(ni - c)n 2 + ni(l — c)] 

-co. - (1 - c) c (32b) 
(n - c)(l - ni) • • • (1 - n fe _i) 
x [-(1 - n k )(n k+ i - c)n k+2 
+ {n k - c)n k+1 (l - c)] 



(1 



k - 1 



(32c) 



In the diagrammatic representation in Eqs. (|32a|l - H32c|l . 
a solid circle • denotes n k+ i or n k+ 2- 

A few words are in order on how to obtain the dia- 
gram of CX given that of X. The Liouville operator £ 
acts much like a differential operator and can be shown 
to follow the product rule C(AB) = A{CB) + {£A)B, 
provided A and B do not involve the same site. Thus, C 
acting on diagrams like those in Table [I] yields a sum of 
terms where C acts on each site individually. Each site 
has one of the diagrammatic elements o, • or - , and 



Cc 



(33) 



6 



element: 




o 


• 


(nothing) 




meaning: 


n 


— c 1 — n 


n 


1 




when averaged, results in: 




1 - c 


c 


1 




horizontal part: 




= 8 


m 


t 


8 2 


meaning: 


(1- 


- n) 2 n(n — c) 


n(l — n) 


n 2 


(n — c) 2 (1 — n)(n — c) 


when averaged, results in: 


(1 


-c) c(l-c) 





c 


c(l - c) -c(l - c) 



TABLE II: Rules for the contributions to a diagram. On top are the elementary diagrammatic elements, at the bottom 
combinations of them. The resulting factors are obtained using Eq. 



Thus, C acting on site i introduces a new diagrammatic 
element • in the diagram on the next site i + 1. If the 
diagram already had an element on that site, one needs 
to multiply the new one with the original one, and this 
gives the three possibilities 



• x o 



= (l-c) 



0. 



(34) 



As the last equation shows, C acting on an element at site 
i yields zero if site i + 1 has the diagrammatic element 
"- " . This is the reason why there are so few diagrams in 
Eqs. (|32c|> : although C could in principle act on all sites 
in Qi(k), most have an element "- " to their right and 
yield zero. 

The matrix elements (Q \£\Q ), (Q \C\Qi(k)) and 
(Qi(k')\C\Qi(k)) are found by taking the inner product 
of the diagrams in Eq. (|32|l with those of Qo and Qi(k') 
in Tabled Because each open circle needs to be covered, 
the only nonzero contributions are for k' = k — 1, k and 
k + 1 , and 



(Qa\C\Q ) = 
(Qi(0)|£|Q ) = 
(Qi(0)|£|Qi(0)) = 



-8- =-c 2 (l-c) 
-88 =-c\l-cf 



(l-c)S? 



(35a) 
(35b) 

(35c) 



(Qi{k)\L\Q x {k)) = -S=^«+(l-c)e==8 

= -c 3 (2-c)(l-c) fc+2 if k >(Q5d) 



and 



(Qi(fc + l)|£|Qi(fc)> 



k 



c a (l-c) fc+3 , (35e) 



where the rules in Table ITTI were used. Using Eqs. (|27|l . 
(|31|l and (|35|l . M E becomes the infinite tridiagonal matrix 



M 1 



-vMl-c) 

Vc(l-c) -1 c^/Y~c 



Cy/T 



-c(2-c) 



(36) 

where diagonal dots denote repetition of the last men- 
tioned expression on that diagonal. 



At this stage, M D (t) in Eq. and M c in Eq. (JT^) 
will be set to zero. The reasons for this are twofold. First, 
it allows an exact solution for spin autocorrelation func- 
tions to be obtained that is in good quantitative agree- 
ment with simulations if the density of up-spins c is not 
too low. Secondly, we will later complete the basis such 
that M D is in fact strictly zero. Eq. I|19l) yields in this 
approximation: 



(37) 



Here, the superscript (1) indicates that the result is only 
a first approximation. Below we will make this into a 
systematic approximation scheme in which further, more 
accurate approximations can be obtained. 

According to Eq. I|21|) . the Laplace transform of the 
correlation C(t) in Eq. (|llf> is the top left element of the 
matrix G, 



C = 



dti 



t C(t) = G 11 . 



(38) 



To perform the matrix inversion in Eq. (|37fl . one uses 
the fact that the inverse of a tridiagonal matrix can be 
performed exactly. In particular, the top left clement 
of the inverse of a symmetric tridiagonal matrix can be 




simulation c=0.8 
single domain c-O.S 

simulation c=0.6 
single domain c=0.6 

simulation e=0.5 

simulation c=0.4 
single domain c=0.4 
simulation c=0.3 
single domain c-O.J 
simulation c=0.2 
single domain c-0.2 
simulation c= 0. 15 
single domain c— 0. 15 



FIG. 1: Results for the single spin time correlation function 
C{t) from the single domain basis {Qo,Qi} [by numerical 
Laplace inversion of Eq. 1144ft in Sec. IV Al using Mathemat- 
ics compared to simulation data (kindly provided by Prof. 
H. C. Anderson). 
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written as a continued fraction: 

_ -l 

at h 
bi a 2 b 2 



52 a3 



a i 



(39) 



a 2 



ft 2 



«3 



Combining Eq. I|3{j|) H39[l , one finds 

z + c \/c(l — c) 



<7« = 



v / c(i - c) z + i -cyT 



-c-\/l — c z + c(2 — c) 



f 



z + c — 



c(l - c) 



n 

^40) 



0+1- 



c 2 (l -c) 



* -r c(2 - c) 2+c(2 _ c) 
The repeating part of this expression is 

c 2 (l-c) 



7 « = 

z + c(2 - i 
This clearly satisfies 



c 2 (l -c) 



(41) 



z + c(2 - c) - 



c 2 (l-c) 



z + c(2-c)-7( 1 ) 
which is solved by 

7W = - {z + c(2 - c) - V4cz+(z-c 2 ) 2 } 



(42) 



(43) 



[Note that the solution of Eq. Q42JI with a plus sign in 
front of the square root in Eq. 143|) does not go as 1/z for 
large z, and is therefore not in agreement with Eq. iJHJ]. 
Inserting this result in Eq. I|4(J|I . one obtains the explicit 
form 



(7« = 



1 



2c(l-c) 



z + 2 - (2 - c)c + ^4cz + (z - c 2 ) 2 

(44) 

The exact correlation function in Eq. (|44|l can 
be Laplace inverted numerically using Stehfest's 
algorithm^ii^ (for example, in Mathematical) . For var- 
ious values of c, the results are shown in Fig. ^ an d com- 
pared with data from simulations on the East model. 
Despite the simple form of in Eq. (|44|) . there is ex- 
cellent agreement between this theoretical result and the 
data for 0.6 < c < 1, reasonable qualitative agreement 
up to c ~ 0.5, while the predicted decay is clearly too 
fast for c < 0.5. 



B. Inclusion of neighboring domains: a complete 
basis 



There is a need to extend the single domain basis be- 
cause it does not properly capture the long time behavior 
of C(t) for c less than 0.5, as seen from Fig. ^ The lack 
of quantitative agreement between theory and simulation 
at long times implies essentially that there is important 
slow behavior in the memory function M D (t) given in 
Eq. Ijl6|l that cannot be neglected. 

This situation is reminiscent of that in fluids. There, 
one starts out describing time correlation functions in 
terms of the linear dependence on the hydrodynamic 
fields of mass, momentum and energy density, i.e., one 
takes these to comprise the set A of Sec. IIII? 8 39 . But 
at lower temperatures or higher densities this does not 
suffice because M D (t) turns out to no longer be a fast 
decaying function. To fix this situation, i.e., to represent 
the missed slow behavior in M D (t), one needs to augment 
the linear basis by vectors orthogonal to it. For this, 
one can take products of A (with proper subtractions to 
assure orthogonality); these additional basis vectors are 
called multi-linear modesM^. The coupling of the linear 
modes to the multi-linear ones "renormalizes" the bare 
values of M D found using only A to M D + Z(z) , where Y. 
is a self-energy. The z dependence of this self-energy is 
such that it can describe slowly decaying behavior such 
as long-time tails. Since the multi-linear modes can be 
interpreted as products of linear hydrodynamics modes, 
this procedure amounts to a nonlinear coupling of hydro- 
dynamic modes and is hence called mode-coupling theory. 

Similarly, if for the East model, the matrix M E is 
taken to be represented at the linear level by Eq. 
where the linear basis (using analogous nomenclature as 
above) composing the set of slow variables is taken to be 
A = {Qo, Qi}, then the memory function corresponds to 
an infinite square matrix represented at the linear-linear 
level that effectively renormalizes the matrix elements of 
M E to M E + M D , according to Eq. (JTIJ). Hence, M D 
takes the role of the self-energy here, and must describe 
contributions to the decay of the spin-spin correlation 
function due tos the projection of the dynamical evolu- 
tion onto a space orthogonal to the linear basis set A. In 
other words, the single domain basis does not span the 
crgodic component and fails to capture all the slow dy- 
namics of the spin fluctuation variable Qo- To represent 
the missed slow evolution, the single domain basis set 
must be expanded to include additional slow modes and 
their coupling to the linear modes must be computed. 

To deduce the appropriate extension of the basis, it is 
helpful to realize, as Fig.^shows, that the decay of C(t) 
that is predicted by the extended linear basis, {Qo,Qi}, 
is too rapid. A reasonable explanation of this particular 
deviation is that in the single domain basis, the final spin 
nk+i in Eq. I|22c|l decays regardless of the spin configu- 
ration to its right. As a result, any slowing down effect 
of a persistent down-spin domain to the right of n^+i is 
missed. It therefore seems natural to try to fix the too 
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?a(0,0) = fl(0) 



QaOi.O) = fl(Ji) 

«1 '2 



Q 2 (Q, l 2 ) = S(Q, h - 1) Q2{h,h) = S(l u h - 1) 

TABLE III: Diagrams of the extension of the domain basis. 
(The notation using R and S is used only in the appendix.) 



rapid decay of the C(t) by augmenting the basis with a 
second down-spin domain, 



1 h 
Qi{h,h) = 7 Tl -j , (no - c) TT (1 - njjn^+i 

h+i2+l 

x II ( 1 - n ^)(^ 1 +i2+2-c), (45) 

j'a=2i+2 

which carries an index doublet of which each mem- 

ber lj can take integer values from zero to infinity, and 

Z 2 {hM)=cV\l-c) 1+ ^ +l ^ 2 . (46) 

As before, unnormalized versions Q2(h,h) — 
Z2(li,l2)Q2(h,h) are defined as well: 

h 

Q2(h,h) = {no - c) JJ (1 - njjn^+i 
j'i=i 

h+h + l 

x [| (l-n j2 )(n il+i2+2 -c), (47) 

j2=h+2 

Diagrammatic representations of the Q2(h, h) are shown 
in Table ITTT1 These are orthogonal to Qo, as well as to Qi, 
since in a diagrammatic representation of (Q0IQ2C1, '2)), 
the trailing open circle (71^+^+2 — c) of Q2 is not covered, 
yielding zero, and in the diagram of (Qi(k)\Q2(h, h)}, it 
is impossible to line up the trailing open circles of Qi(fc) 
and Q2(h,h) without having the solid dot • (=n) over- 
lap with a horizontal line - (=1 — n), which yields zero 
[cf. Table ITTj . It is easy to establish that the Q2(h,l2) 
are also orthonormal among themselves. 

There is no obvious reason to stop this procedure at the 
two-domain, or "bi-linear" , level and, in fact, the basis 
can be extended to a complete set in the relevant ergodic 
component in a straightforward fashion. The elements of 
this complete basis are most easily written diagrammat- 
ically as a sequence of a down-spin domains of different 
sizes kj, separated by single up-spins: 



t (&1 : ■ ■ ■ 1 k a ) 



1 



Z a (ki , ■ ■ ■ , k a ) 



Here, a = . . . 00, kj = Q . . . 00 (j = 1 . . .a) and 



(48) 



Z a {ki 



1 k a ) 



c i/2 (1 _ c) i/2 ifa = 

c (i+ a )/2^ _ C )i+E?=i fcj/a otherwise. 

(49) 



It is easy to see that the Q a are all independent: the inner 
product of two of them is zero unless their diagrams [as 
in Eq. jlBjl] are equally long, so that both open circles are 
covered. But then the interior of the diagrams has to have 
matching top and bottom parts too, otherwise a solid 
dot and a horizontal line occur at the same site, and this 
gives zero (Table Hljl . The only nonzero inner product of 
a Q a (ki, . . . k a ) is therefore with itself. Due to our choice 
of normalization ([Q a (ki, . . . , k a )] 2 ) — 1. Since each Q a 
is orthonormal to all others, each contributes a unique 
direction in the Hilbert space to the basis which could 
not be formed from the others: The Q a are independent. 

To also establish completeness, we will now count the 
number of element of the above basis in the finite system 
of spins. In that system, there exist only Q a for which 
a+ 1 +X^=i kj < N (which also limits a < N). Elemen- 
tary combinatorics shows that the number of different Q a 
for given a and N is ( N ^ 1 )- The total number of Q a is 

thus ~ X ( N ~ X ) = 2 JV " 1 . Thus the above set of 
basis vectors covers only half of the full Hilbert space, 
which has 2 N dimensions. But it is easy to see which ba- 
sis vectors are missing and why they are not important. 
The expression in Eq. 148(1 always starts with uq — c, even 
though the first spin can have two values. Independent 
vectors can be found by taking a different expression for 
the first spin. In fact, one can take 1, i.e., one could con- 
sider variants of the basis vectors in Eq. I(48|l in which the 
factor of (no — c) is not present. Call these Qp, of which 
there are as many as there are Q a . These new vectors are 
orthogonal to each other as well as to the Q a because in 
{QaQp) the initial o in the diagram of Q a is not covered 
by the diagram of Qp and (o) = 0. Thus, the Qp are 
the missing basis vectors. However, they are completely 
unimportant here because the Qp are also orthogonal to 
HQci-, as is seen from the fact that CQ a will always have 
a "o" as a first element (£0 = —cm ) so that the inner 
product with Qp is zero (for this it is important not to 
have periodic boundary conditions). 

So the basis set Q a is not a complete basis for all possi- 
ble spin configurations, but it is a complete orthonormal 
basis for all spin configurations to which the Q a couple. 

These considerations also imply that the East model is 
not ergodic: the state space contains at least two ergodic 
components, which are such that a configuration in one 
of them can never make a transition to any configura- 
tion in the other. Noting that the space spanned by Qp 
contains all quantities insensitive to the value of no, one 
realizes that it constitutes an East model with an effec- 
tive length of A^ — 1. The argument above then shows 
that the state space of this smaller East model can also 
be split into at least two ergodic components. Applying 
this argument recursively reveals that there are N + 1 er- 
godic components. The p-th ergodic component consists 
of functions not sensitive to the values of spins no through 
some rip-x, with < p < N, and has 2 N ~ P ~ 1 dimensions 
if p < N and one dimension if p = N. The collection of 
these ergodic components has 1 + Y^p^Q 1 2 N ~ P ~ 1 = 2 N 
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dimensions, and thus indeed spans the full Hilbert space 
of the spin chain of length N. 

Since we are interested in the time auto-correlation 
function of spin no, the relevant ergodic component is 
the one spanned by Q a , and we conclude that the Q a 
are the only basis vectors needed. Having established the 
"relevant completeness" , one can take the limit TV — > oo 
again, so we need not worry about the boundary condi- 
tion imposed on npf. 

The extension of the basis set to include an arbitrary 
number of domains is useful in developing a systematic 
approach to generate successive improvements for C for 
lower c. Since the basis A = {Q a } spans the ergodic 
component of Qo, it follows that ip(t) = e( 1 ~ v * >ct (l — 
V)LQo = and the memory function M D (t) vanishes. 
Hence, from Eq. (JTHJ), G(z) = (zl - M)~\ where the full 
matrix M = M E in this complete basis can be written as 



M 00 M i 
M 01 Mn M 12 

M\ 2 M 22 



(50) 



where M a p = {Q a \C\Qp) and it was used that this is zero 
unless \a — (3\ < 2, as can easily be shown (also, diago- 
nal dots do not denote repetition now: all M Q/ 3 can be 
different). Generalizing Eq. by repeatedly applying 
the matrix equality (see e.g. Ref— , page 70) 



a c 
d b 



one finds 



:dl-l 



-[a - -^-icb- 1 



C(z) 



[b-^da- 1 [b-^] 



zl-Moo -M i 
-M i zl-Mn -M12 

-M} 2 zl-M 2 



dcl-l 



(51) 



(52) 



li 



zl 



1oi Mj x 



-(53) 



100 



zl 



12 



Zl - 



zl 



100 



1oiM 01 



(54) 



zl-Mn-En(z) 



where the self-energy matrix at the linear-linear level is 



defined to be 



Mi 9 M 



t 



12 'VI 12 



zl - M 22 



M 23 Mj 2 



zl - M? 



M12 M\ 2 



zl - M 22 - S 22 (z) 



(55) 



(56) 



For convenience a non-standard (but unique) notation 
for a matrix fraction has been introduced here, such that 
if A, B and C are matrices then 



AB 
~~C~ 



A-C 



B. 



(57) 



This notation saves a lot of space and avoids many nested 
parentheses and inverses that would be required in more 
standard notation. We remark that Eq. (|52|) is similar to 
the matrix formalism of Wu and Caoi^, while Eq. (|53|) 
has similarities to the continued fraction formalisms of 
Mori 47 i 48 and Schneider—. The structure of Eq. lj52*|) is 
that of a mode-coupling theory in which the role of mode 
order is played by the number of domains. The effect of 
the higher order modes is to renormalize the "transport" 
coefficients approximated by M E at the linear level. 

By truncating Eq. I|53|l at ever deeper levels, i.e., set- 
ting M a a+ i = 0, corresponding to T, aa (z) — 0, for in- 
creasing a one gets expressions which work well for ever 
lower values of c, denoted as & a '{z). For example, trun- 
cating at the zeroth level gives 



C<°\z) 



1 



1 



M 



00 



I 1 

z + c 



(58) 



while truncating at the first, linear, level yields the result 
in Eq. I|44|) . Following this procedure further, the first 
correction to the linear basis results involves evaluating 
the self-energy in the approximation where one ignores 
the effects of three-domains and higher, i.e. [cf. Eq (|56|) ]. 



Eii(z) 



Mi 2 M 

zl 



'122 



(59) 



corresponding to a bi-linear type of mode-coupling the- 
ory. Due to the simplicity of the coupling with respect 
to mode order and for different domain sizes, an exact 
expression for this approximate self-energy can be ob- 
tained. In the appendix £n(z) is explicitly evaluated to 
be: 



En 



-VI - cfji 



VI - cfji (l-c)?7i+772 ~ cm 

-VI - crj 2 (2-c)fj 2 



(60) 
and 



where the functions fjj are given by Eqs. ljA.21 
(|A.25I> . Using this expression for the self-energy, the 
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linear-linear matrix G(z) of the previous section (which right hand side of Eq. I|52|) incorporating the zeroth and 
is in fact the top-left block of the inverse matrix on the first level) is renormalized to 



J 




z + c a/c(1 - c) 

a/c(1 - c) 2+I-771. 



-VI - c(c- 771) 



VI — c(c - 77i) 2:+ (2 - c)(c — 772) + (1 - c)(772 - fji) -y/l - c(c - 7)2) 

-VI - c(c- 772) z + (2 - c)(c - 772) 



(61) 



(62) 



from which the single spin time correlation function = [G^ ]n is computed with the continued fraction expression 
in Eq. ^ to be 



2) 



1 



where a(c, z) is 



a(c, z) = z + 1 — 771 



a(c, z) 
{l-c){c-f,xf 



(63) 



In Eq. I|64|). the repetitive part 7^) satisfies 



2+ (2 - c)(c - 772) + (1 - c)(772 - 771) - 7(2)' 



- (2) = (l- C )( C -^ 2 ) 2 



z + (2 -c)(c- 772) - 7 (2) ' 
This is solved by 

7 (2) = ^ + (2 - c)(c - 772) - ^4(c - 772)2 + (z - c(c - t^.)) 2 } 
Note the resemblance with 7' 1 ) in Eq. I|43|l . 



(64) 



(65) 



(66) 



As before, the exact result in Eq. (|63|1 is inverted nu- 
merically using Stehfest's algorithm^MI^. For various 
values of c, the results are shown in Fig. [3 and compared 
with data from simulations of the East model. Notice 
that there is a huge improvement over the results ob- 
tained using only the single domain basis {Qo,Qi} i n 
Fig. ^ There is now excellent agreement between the 
theory and the data for c > 0.4, and reasonable agree- 
ment up to c sa 0.3. Furthermore, while the theoretical 
decay is still too fast for c < 0.3, the small time behavior 
is captured beautifully. In particular, the shoulder that 
appears in the simulations for low c is reproduced by the 
extended theory as well, something the single domain ba- 
sis could not do. 

In the low temperature region (small c), the long 
time behavior of the spin autocorrelation function pre- 



dicted by the two-domain basis set is well-described by a 
stretched exponential C(t) ~ exp [— {t/r) 13 ] with a tem- 
perature independent stretching exponent of /? w 0.6. Al- 
though it is encouraging that the stretched exponential 
time profile is indeed predicted by the theory, simulations 
indicate that in fact the stretching exponent j3 should 
have a weak temperature dependence^, with j3 decreas- 
ing in value as the temperature decreases. The origin of 
this discrepancy between our theory and numerical sim- 
ulation is not clear and is under investigation. 

In principle, the effect of three down-spin domains (tri- 
linear modes) can be included in the same spirit, i.e., by 
evaluating the self-energy at the two-domain level 222(2) 
using matrix methods similar to those applied to obtain 
Eq. (|60|) . Unfortunately, the algebra becomes even more 
cumbersome and explicit evaluation of the self-energy 
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simulation c=0.8 
extended basis c=0.8 

simulation c=0.6 
extended basis c=0.6 

simulation c=0.5 
extended basis c=O.S 

simulation c=0.4 
extended basis c=0.4 

simulation c=0.3 
extended basis c=03 

simulation c=0.2 
extended basis c=0.2 
simulation c=0.15 
extended basis c=0.15 



uous, so the more general expression to replace Eq. (|67() 



is 



10000 100000 



FIG. 2: Results for the single spin time correlation function 
C(t) using the extended basis {Qo,Qi,Q2} in Eq. 1631 of 
Sec. IV Bl (numerically Laplace inverted using Mathematica), 
compared to simulation data. 



matrices at higher and higher order becomes effectively 
impossible. Alternatively, one can resort to numerical 
approaches in which the maximum domain size k m is 
fixed and all matrix inversions are carried out numeri- 
cally. By monitoring convergence to a set level of preci- 
sion, such a procedure provides a systematic and numeri- 
cally tractable method of predicting the decay of the spin 
autocorrelation function for arbitrary values of c. 



VI. RELAXATION BEHAVIOR 

One of the main advantages of the matrix method out- 
lined here is that it is straightforward to obtain analytic 
predictions for rather detailed features of the dynamics. 
For example, one of the commonly calculated quanti- 
ties from simulation data is the relaxation time r. For 
systems exhibiting such non-trivial relaxation behavior 
as stretched-exponential, the definition of the relaxation 
time is a matter of choice. Perhaps the most sensible 
way to view the relaxation time for such systems is to 
consider it as the weighted- aver age of a distribution of 
relaxation times. For example, based on the spectral de- 
composition of the Liouville operator, one can formally 
write the spin autocorrelation function as a weighted sum 
of exponentials with relaxation times r„ , 



C(t) = 2J c n exp(-i/ T„ 



(67) 



Since the Liouville operator is Hermitian and the spin 
variables are real, one is guaranteed that the relaxation 
times r n and coefficients c n = (Qol^n) 2 , where \ip n ) are 
the right (and left) eigenvectors of the Liouvillian £, are 
real and positive. Furthermore, since C{t = 0) = 1 = 
XL c nj the coefficients c n are proper weights for the re- 
laxation time r n . However, since C is of infinite dimen- 
sion, its spectrum can be (partially or completely) contin- 



C(t)= J P (T')exp(~t/r')dT\ 



where p(r') > 0, p(r' < 0) = and / p{r') dr' = 1. One 
can therefore define the average relaxation time as 



t = J P (t')t' dr'. 



(69) 



Noting that the Laplace transform C(z) of Eq. I|68|l is 

J Z + l/T 1 

we see that 

t= / p(t')t' dr' = C{z = Q). (71) 



Note that in the case in which a single relaxation time 
t* dominates all others, one observes that r w r* since 
p(r')^6(r'-r*). 

Note also that in taking the point z = 0, the expression 
is sensitive to long time behavior. This in contrast to e.g. 
the average rate J p(r')(l/r')(iT / which by Eq. (|68|) is just 
—(d/dt)C(t = 0) = c and contains no information on the 
long time behavior. 

Given the analytical results for the Laplace transform 
of the spin autocorrelation function in the one-domain 
[Eq. (|44|l ] and two-domain representations [Eq. ^63(1 ] of 
the slow dynamics, explicit expressions for t(c) can be 
obtained by setting z — in the respective equations. 
For example, not including domains as in Eq. (|58|l gives 
r™ = 1/c, while in the one-domain basis, Eqs. I|44|) and 
(|71|l lead to the simple result 



T« = 



1 - c- 



(72) 



in which the average relaxation time diverges as c~ 3 as 
the concentration c approaches zero. Furthermore, in the 
two-domain representation, the average relaxation time 
is a complicated function of c. In the limit that c — * 0, 
we find that ~ c~ 4 . In Fig. ©, the theoretical 
predictions of the average relaxation time in the one- 
domain and two-domain basis sets are compared with 
numerically-integrated simulation data. Note that as is 
evident from Figs. QJ and the two-domain predic- 
tions significantly improve the one-domain results but 
still underestimate the relaxation time of the system at 
small values of c. 

From the relationship between c and (3p in Eq. J2J), it 
is clear that at low temperatures c ~ exp(— (3p). Since 
the logarithm of the average relaxation time log(r) is 
proportional to log(l/c) for t^°\ and t^ 2 \ a plot 
of log(r) versus (3p, yields a straight line in the small c 
(low temperature) limit. Thus we can conclude that the 
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FIG. 3: Logarithm (base 10) of average relaxation time r for 
various values of c. The inset shows the same as a function 
of the logarithm of c (also base 10) for the theoretical results, 
with slopes of — 1, —3 and —4 showing their scaling behavior. 



zero, one and two-domain basis sets all yield a relaxation 
time that diverges according to the Vogel-Fulcher law 
t ~ exp{— const/(T — To)} with a glass transition tem- 
perature of T = 0. Note that these results are in contrast 
with the exact result for the equilibration time r e of a sys- 
tem quenched to T = where r e ~ exp{const/T 2 }^2*iS. 
This finding is somewhat surprising given that the equi- 
libration time was calculated in the asymptotic small c 
regime using ideas of domain structure rather similar to 
those presented here. 

Given the relatively simple structure of the matrix 
G(z), it is easy to numerically examine many detailed 
features of the relaxation given a finite domain basis set 
specified by setting a maximum domain size k m . For ex- 
ample, one can easily examine how the spectrum of C 
depends on c. At the same time, the actual distribu- 
tion of the c„ can be computed numerically to see how 
many relaxation modes are relevant as a function of c. 
From this information, one can try to attempt to estab- 
lish a link between the distribution of relaxation times as 
a function of temperature and the asymptotic stretched 
exponential form, as suggested in reference^. This may 
be an instructive way to examine the failure of the two- 
domain basis to correctly predict the temperature de- 
pendence of the stretching exponent [3. However since 
analytical results are available for all quantities, it is de- 
sirable to obtain analytical expressions for such features 
as the width or spread a of the relaxation times r' as a 
function of c. The spread in r' is defined by 



p(T')(T' - T) 2 dT'. 



(73) 



Now one can use that C"(0) = Um*_»o d/dzC{z) — 
— J p(t')t' 2 dr 1 [cf. Eq. (□}] to write a = 
{-C"(0) - [C(0)} 2 } 1/2 . Since we have obtained closed 
expressions for C(z), analytic expressions can be ob- 



tained for a. For example, using the one-domain basis 
set, we find that 



(74) 



whereas the expression for in the two-domain basis 
is a complicated function of c [note: c' ' is actually zero]. 
From these analytical expressions for a, one immediately 
sees that, in fact, a diverges as c approaches zero in the 
same way as r. Hence a plot of cr/r remains finite for all 
values of c. Furthermore, noting that in the one-domain 
basis, 



l-c + c 2 



(75) 



it is evident that lim c _»o a^- 1 ' /t^- 1 ' = 1. Surprisingly, the 
same conclusion holds in the two-domain basis, as is ev- 
ident from Fig. Note that at large values of c w 1, 
a ~ indicating that the relaxation is dominated by a 
single mode. 

We note that higher order derivatives of C(z) at z = 
can similarly be used to investigate further characteristics 
of the relaxation time distribution such as the skewness 
and the kurtosis. More generally, Eq. I|tj8fl shows C(t) to 
be the Laplace transform of the distribution of relaxation 
rates. If r = 1/r are the relaxation rates, then their 



distribution is P(r) 
written as 

C(t) = 



2 p{l/r) and Eq. JfjEJ) can be 



P(r) exp(— rt) dr. 



(76) 



In this sense, C(i) is the Laplace transform of P(r). 
Thus, given C(t), one might expect to be able to use 
the numerical Laplace inverse of the Stehfest algorithm 
to obtain P(r). Unfortunately it turns out that using 
Stehfest's numerical Laplace inverse method on C(t), 
which was itself obtained from C{z) by the same method, 
is unstable; it in fact yields an incorrect result for P(r), 
namely a highly oscillating function, which is not non- 
negative and not normalized to one. Since the distribu- 
tion of relaxation rates was considered here mainly as 
an illustration of the power of the theoretical approach 
presented in this paper, solving the numerical instabil- 
ity associated with applying the Stehfest algorithm twice 
and determining the distribution of relaxation rates in 
detail, is left for future work. 



VII. HIGHER ORDER CORRELATION 
FUNCTIONS 

Given the glassy nature of the dynamics of the East 
model, it is interesting to probe higher order correlation 
functions to examine issues of cooperativity in the dy- 
namics and non-Gaussian statistics. In particular, one 
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FIG. 4: Ratio cr/r of the width a of relaxation time distribu- 
tion to the average relaxation time r as a function of c, in the 
one-domain and two-domain bases. 



can look at the neighbor-pair spin correlation function 

(n i (t)n i+ i(t)n i (0)n i+ i(0)) = (Qi(fc; *)&(&; 0))4,o 

= G 22 (t), (77) 

and a related quantity 

A(i) = G 22 (t) - {hi(t)hi(0)){h i+1 (t)h i+1 (0)) (78) 

that examines the non-Gaussian nature of the normalized 
spin fluctuation variable hi. Given the simplicity of the 
matrix method, it is relatively straightforward to obtain 
analytic expressions for higher order correlation functions 
such as Eq. I|77|) . For example, from the definition of the 
neighbor-pair spin variable, which corresponds to the lin- 
ear basis set element Qi(0), it follows that the Laplace 
transform G 22 of the function G 22 (t) is the 2 — 2 element 
of the infinite matrix G, which, in the two-domain ba- 
sis approximation, is given by Eq. (|61|l . Using standard 
matrix inversion methods, the 2 — 2 clement of is 



r<( 2 ) _ 

U 22 — 



1 



a[c, z ) 



e(l-e) ' 
z-\-c 



(79) 



where a(c, z) is given in Eq. i|tj4|) . In Fig. JSJ, the func- 
tions G 22 (t) and A(t) are plotted versus time for vari- 
ous values of c (using Stehfest's algorithm for the inverse 
Laplace transform). Note that the agreement between 
the theoretical predictions and the simulation data is ex- 
cellent for all times for all but the smallest value c = 0.2. 

The neighbor-pair autocorrelation function exhibits 
several remarkable properties that are rather unlike those 
of the spin autocorrelation function. Note that in the 
short time limit t < 1 the relaxation of G 22 (t) is in- 
dependent of the equilibrium up-spin concentration c. 
This result can be explained by examining a short time 
expansion (large z) of G 22 , from which it is seen that 



a(c, z) ~ z + 1 and hence G 22 ~ 1)' correspond- 

ing to simple exponential relaxation G 22 (t) w exp(— t). 
Effectively this approximation corresponds to the short 
time expansion G 22 ~ l/(z — (Qi(0)|£|Qi (0))). Even 
more remarkable is the clear emergence of a plateau in 
the neighbor-pair autocorrelation function as c decreases 
and the system becomes "glassy" , yielding a two-step 
relaxation time profile similar to that observed for the 
dynamic structure factor at microscopic length scales in 
simple glass-forming systems. In such systems, the on- 
set of the plateau, generally called the /3-regime, is rela- 
tively insensitive to temperature and is often associated 
with the phenomenon of dynamic caging in dense fluid 
systems. In this regime, fluid particles typically oscil- 
late in the traps formed by their immediate neighbors 
and little relaxation of the system occurs. This behav- 
ior typically continues until a typical time scale, known 
as the a-regime, is reached in which particle cages are 
temporarily broken. This a time scale is strongly tem- 
perature dependent and scales with the overall relaxation 
time of the system. Interestingly, similar behavior is ob- 

served in G 22 (t) of the East model: There is an initially 
rapid decay (with time scale t ~ 1) at which point a 
plateau appears. The plateau typically extends to times 
corresponding to the average relaxation time r of the spin 
autocorrelation function. However, unlike simple liquid 
systems, the height of the plateau is strongly temperature 
dependent, occurring roughly at value of c. In the East 
model, one can interpret the emergence of the plateau as 
arising from a kind of effective dynamic caging of the pair 
spin variable n^ni+i that occurs when ni + \ = 1. When 
the right neighbor of a given spin i is up, the spin 
can oscillate between values of 1 and for extended pe- 
riods of time, corresponding to a kind of vibration in a 
cage. This behavior will persist until the spin i + 1 flips, 
which typically will occur at times t ~ r. Furthermore, 
the probability of finding such a caged system scales with 
the likelihood of finding an up-spin in equilibrium, c. 

The two-step relaxation of G 22 (t) was also found nu- 
merically by Wu and CaoA^ (who refer to this quantity 
C2). Wu and Cao showed that the relaxation can be 
described with a stretched exponential behavior at long 
times. From the numerical analysis of our theoretical 
expressions, we find that the parameter-free theoretical 
relaxation profile is also well-described by a stretched ex- 
ponential with the same stretching exponential (3 ps 0.6 
found in the analysis of the spin autocorrelation function 
C(t). 

As can be seen from Fig. (JSJ), the spin fluctuations hi 
do not behave as Gaussian random variables at all time 
scales and for all values of c, unlike their counterpart, 
the Fourier components of the mass density, in simple 
liquid systems. It can also be observed that the decay 
of the spin fluctuations is slower than that predicted for 
a system exhibiting Gaussian statistics for all times at 
high values of c. As c drops below 0.5, the decay be- 
comes faster than Gaussian at short times but slower 
than Gaussian at long times. The fact that c = 0.5 is 
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FIG. 5: Higher order correlation functions. Top: The 
non-Gaussian measure A(t) defined in Eq. 1781 . Bottom: 
The neighbor-pair auto-correlation function G22(i). In both 
graphs, open circles and solid lines correspond to theoretical 
and simulation values, respectively. 



special can be seen from a short time expansion of A(t): 

A(t) = (&(0)|e £ *|&(0)) - (Q \e ct \Qo) 2 

= 1 + *<&(0)|£|&(Q)) + ^-(QxWI^IQiCO)) 



1 + i(Q |-C 2 |g > + ^-(Q |/: 2 IQo> 



0(t 3 



Using the rules elaborated in Sec. IV Al all quantities ap- 
pearing above are easily evaluated to reveal the exact 
result: 



A(t) = (2c- l)t 



l-(2c+l)| 



+ 0(t 3 ), (80) 



from which the sign change for c = 0.5 is explicitly evi- 
dent at short times. 

One can also note in Fig. (jSJ) that the maximum pos- 
itive deviation from Gaussian behavior (i.e. slower than 
Gaussian) occurs at a time which scales roughly with the 
average relaxation time r. 



VIII. DISCUSSION 

In this paper, the East model — a linear kinetically 
constrained spin model which is statically structureless 
— was studied theoretically taking the domains of down- 
spins as a starting point. The constraints in the model 
lead to a very slow spin relaxation for low up-spin density 
c because of the existence of these down-spin domains, 
of which each spin has to flip at least once before a spin 
on the left of the domain can relax. Such highly coop- 
erative, hierarchical events driving the relaxation mimic 
heterogeneous behavior in glasses. 



The way the down-spin domains were taken into ac- 
count was by using them in the construction of a basis 
which is complete on the relevant ergodic component. In 
the complete domain basis, the theory is formally exact, 
but the basis needs to be truncated to get explicit results. 
In this truncation, one only limits how many simultane- 
ous domains are included without restricting the possible 
sizes of those domains. When we restricted ourselves to 
a single domain description, an exact result for the single 
spin time correlation function C(t) [C^] was obtained 
which gives a good quantitative description for c larger 
than about 0.5. An extension including neighboring do- 
mains led to an exact expression [C^] which described 
the slow, glassy behavior correctly down to c « 0.3. A 
general procedure was outlined to obtain further approx- 
imations. 

The main advantages of our approach over others are 
that a) it gives explicit analytical results without fitting 
parameters, b) it requires neither an arbitrary closure for 
the memory kernel nor the construction of an irreducible 
memory kernel such as in mode coupling theories and c) 
nonetheless, it described low c behavior equally well as 
these mode coupling theories. The explanation for this 
power is that domains of all sizes are included. 

At a given level of truncation, the matrix approach out- 
lined here allows analytical results for the spin autocorre- 
lation function to be obtained. Armed with these results, 
it is possible to assess the effect of truncation the multi- 
domain basis by evaluating approximate expressions for 
the "self-energy" terms, as was done in Sec. [V] One can 
then examine the time scale at which the higher domain 
corrections become important and their magnitude for a 
given value of c. Such information is useful in examining 
dynamical scaling relationsiS^i. 

The matrix approach is also well-suited for examining 
higher-order correlation functions, such as the neighbor- 
pair auto-correlation function, that probe detailed as- 
pects of the dynamics, as was shown in section IVIII 

Our theory does not require an ansatz for a closure 
relation between the memory kernel and the correlation 
function, yet it does have the structure of a mode cou- 
pling theory. First of all, the theory, derived using a 
projection operator formalism, yields a basis set very 
similar to the multi-linear set in the theory of Oppen- 
heim et al. Secondly, successive truncations of the set are 
like including only linear modes, or also bilinear modes, 
or also tri-linear ones, etc., again very similar to mode 
coupling theories for fluids. Finally and perhaps most 
strikingly, without assuming a closure relation, a self- 
consistent equation emerges for part of the result, i.e., 
for 7« in Eq. of Sec.|Vl]and for in Eq. (JHSJ 
of Sec. lVBl Thus, in a sense, the correct closure relation 
follows unambiguously from the theory rather than being 
assumed. Perhaps this is an indication why mode cou- 
pling theories can work, at least in some range of c, if the 
closure relation is well chosen. However, as the difference 
between the closure for 7W and "j^ 2 ' shows, the required 
closure depends on how low c is. The closure can also be- 
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come "hierarchical" , in the sense that depends on 772, 
which itself satisfies a self-consistent equation. 

A natural question is how adaptable is the matrix ap- 
proach outlined here for other conditions of spin facilita- 
tion, such as the Frederickson-Andersen&i (FA) model, 
higher-dimensions and other types of lattices. The exten- 
sion to the FA model involves extending the one-domain 
basis set to include domains on both sides of the targeted 
spin and involves slightly more complicated matrix alge- 
bra than that presented here^. For the FA model for 
which self-consistent closure schemes in the context of 
mode-coupling theory appear to work quite wellS, quite 
good quantitative agreement can be obtained with the 
simple single-domain basis set. Extensions to include 
multiple domains can be carried out numerically for the 
FA model as well as other generic models. In addition, 
higher-dimensional models can also be tackled in a nu- 
merical fashion using finite basis set representations, pro- 
vided the basis sets include domains that are sufficiently 
large. Although finite matrix representations are always 
bound to give the incorrect long time asymptotic behav- 
ior for systems exhibiting stretched-exponential profiles, 
the short and intermediate time behavior can be repro- 
duced with great accuracy. 

It is conceivable that the complete basis set presented 
in Sec. IV Bl has a deeper structure that could be exploited 
for the description for c — > 0. Also, the domain basis 
might be used to describe the response of the East model 
to a sudden "quench" to low c values. Work on these 
issues is in progress. 

Finally, our approach shows how important it is to 
first identify the "slow modes" of a system, in this case 
the down-spin domains, before embarking on a mode 
coupling- like description of the long time behavior of cor- 
relation functions. 



it is clear that we must evaluate matrices such as M12 = 
(Qi(l)\£\Q2{h,h))- The double indices on Q2{h,h) 
tend to make the algebra somewhat less transparent than 
in Sec. IV Al and it turns out that the self-energy matrix 
can be evaluated more easily by splitting up the set Q2 
into two-domain variables for which I2 — and those for 
which I2 > 0, by defining 



R{k) = Q 2 (fc,0) 

s(h,h) = &0i,i 2 + i) 



(A.l) 
(A.2) 



and likewise for unnormalized versions, as indicated in 
Table imi The matrix M 12 then takes on the form M 12 = 
[M QR , M qs ], where 



Mqr = (QAC\R) 
M QS = (Qi\£\S) 

and M22 is written in the block form 



M 



22 



M] 



where Mrr, 



•RS ™SS 

Irs and Mss are 

Mrr = (R\C\R) 

Mrs = (R\C\S) 
M ss = (S\C\§), 



(A.3) 
(A.4) 



(A.5) 



(A.6) 
(A.7) 
(A.8) 



and where the notation that R or S without any argu- 
ment denotes the column or row vector composed of all 
R(k) or S(ki,k2) respectively. 

The matrix Mqs is in fact zero, so the matrix self- 
energy En is 
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En 



M 



QR 



Zl -M 22 ] RR M QR . 



(A.9) 



Using the matrix equality Q5ip. this expression can be 
re-written as 



En = M 



QR 



Zl - Mrr 



Mqr M^ fl 



M rs M 



RS IVI flS 



zl - M ss 



Mqr 



zl - Mrr 



Mrs M rs 



(A.10) 



zl 



APPENDIX: SELF-ENERGY MATRIX IN THE 
TWO-DOMAIN BASIS 

From the expression for the self-energy matrix in the 
two-domain approximation, 



En(z) 



M12 M\ 2 
zl - M 22 ' 



The explicit calculation of all the matrix elements ap- 
pearing in Eq. (jA.lOfl proceeds as follows: We start with 
Mqr defined in Eq. (|A.3|I . Combining the diagrams of 
CQi(k) in Eqs. I)32b|l and H32c[l with the diagrams of 
R(k') = Q 2 (fc',0) in Table |III| yields 

(Q 1 (k)\£\R(k-l)) = (l-c)S =c 3 (l-c) fc+3 
(Qi(k)\C\R(k)) = -« =-c 3 (l-c) fc + 3 



1G 



while all other (Qi(k)\C\R(k')) are zero. By similar dia- 
grammatic means, one finds (S(h, l 2 )\L\Qx{k)) =0 — so 
that Mqs in Eq. (IA.4|) is indeed zero as we anticipated 
above. Also (R(k)\£\Q ) = (S(h,l 2 )\£\Q ) = 0, con- 
firming that M 2 = 0. Using Eqs. (|A.3jl . as well as JH|) 
and gives 



M 



QR 



'-(l-c)c 1 / 2 
(l_ c )3/2 c l/2 



(A.ll) 



Next we will determine Mrr defined in l|A.6|) . For this 
we need the diagrams of £R(k): 



£R(0) = -c 
CR{k > 1) = -c 



— (1 — c)ccm — 




(A.13) 



Combining these diagrams with those of R(k') in Ta- 
ble HTH yields 



(R(p)\C\R(0)) = -8B» -(l-c)8B -8B 

= -c 3 (l-c) 2 (2-c + c 2 ) (A.14) 

(R(k)\C\R(k)) = -(l-c)8=88 + S=±S3 

= -c 3 (l-c) k+2 (l + c 2 ) if k> 1 (A.15) 

while (R(k)\C\R(k')) = if k £ k'. Using also Eqs. gBJ 
and l|47[) . one finds 

[ M RR\kk' = -( 2 - c + c 2 )4o4fc' - (1 + c 2 )(l - 4o)4fc', 

Next to determine is M rs ■ Combining the diagrams 
of S i n Table EH with those of CR in Eqs. (|A~12)l and 
(|A.13fl . we see that many elements of Mrs are zero, while 
the nonzero ones are restricted to 

(S(k, 0)\C\R(k)) = -&4n2S = c 4 (l - c) k+3 (A.17) 
for all k > 0. Using Eqs. gUJl and (|A~7)l . we find 

[M RS } kllh = c(l - c) 1 ' 2 8 hk 8 h0 (A.18) 

The final matrix to determine is Mss, for which we 
require CS(h, l 2 ): 



CS(0,l 2 ) 
CS(h > l,h) 



h + 1 , n , z 2 

-c» o + (1 — cjc« 0» 



h + 1 

c« 0» 



h-1 h + l s h h 

c □» o + (1 — cjc ■ 0» 

h h + 1 
— o ■ c» . 



while all other (S'|£|5) are zero, Combining with 
Eqs. (gOJl and (|A.8|I . one gets 

[MasW^j = V[(^i 2 +i + ^ 2 -i)c(l-c) 1/2 
-^i 2( 5; l0 (l + 2c-c 2 ) 
-*jji,(1-*Jio)c(3-c)]. (A.19) 



In view of Eqs. (|A.10|I and (|A.18fl . we need the l 2 = 
0, 1' 2 = component of the inverse of zl — M55. This 
matrix is diagonal in l\ and l[ and tri-diagonal in l 2 and 
1' 2 for fixed l\ and l[. Thus, we can use Eq. i|39|) to write 



M 



friO 

Si - £i 



^, (A.20) 



«2 - £2 



where 



a 1 

«2 



z + 1 + 2c - c 2 
z + c(3-c). 



(A.21a) 
(A.21b) 



and £i and e 2 result from the repeating part of the contin- 
ued fraction that results from applying Eq. I|39l) . Similar 
to 7W in Eq. (|4~2)l in section IV^l they satisfy 



c 2 (l - c)/{a,j - ij). 



(A.22) 



With the requirement that they go as 1/z for large z, the 
solutions are 



,-.^ 2 -4 C 2 (l- C 
- 1 M t 



(A.23) 



The subexpression Mi?s[zl — M$s 
now becomes, using Eq. (lA~18ll . 



in Eq. ||A~T0|) 
20j, (|A3T|) and 



Mrs M\ 



RS 



Zl 



<ss 



= Skk> [ei4o + £2(1 - 5, 



k(>) 



(A.24) 



J kk' 



Since this matrix and M55 in Eq. I|A.16|) are diagonal, 
the inverse of (zl 
simply 



M RR - M RS {zl - MssY-^rs) is 



a- M RR - 



M R sM rs 



-1 -1 



zl- 

SkO^kk' 



<SS 



J kk' 



(1 - 6*0)6. 



kk' 



2 - c + c 2 - ei 



Skk' 



■>ko 



1 + c 2 - e 2 

1 - 4o 

2 /5T„ 



l-2c + c 2 /£i l-2c + c 2 /£ 2 



Combining with the diagrams of S in Table ITTll one finds 

(S(0,l 2 )\C\S(0,l 2 )) = -(l + 2c-c 2 )c 3 (l-c) 3+h - 

(S(0,l 2 + l)\C\S(0,l 2 )) = c 4 (l- c ) i2 + 4 

(S(h>lJ 2 )\C\S(hJ 2 )) = -(3-c)c\l-c) h+l * +3 

(S(h>l,l 2 + l)\C\S(hJ 2 )) = c 4 (l-c)^+ 4 , 



where in the last equality we used again Eq. l|A~22)l . 
Given this last form, we can shorten many equations by 
using the expressions ijj = c(l — c)/(l — 2c + c 2 /ej), which 
are explicitly given by 



m = 



c(l - c) 



1 - 2c 



2c 2 



(A.25) 



4c 2 (l-c) 
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and in terms of which we have 



zl 



'RR 



Mas_M 
zl - M 



TVS 

ss 



kk> 



Skk'lmfiko + 772(1 - 4o)] 

c(l - c) 2 



Inserting this result in the expression for the self-energy 
matrix in Eq. (|A.10I) , one obtains the result presented in 
Eq. 
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